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ABSTRACT 

This paper presents an initial survey of the properties of accretion flows in the 
Kerr metric from three-dimensional, general relativistic magnetohydrodynamic simula- 
tions of accretion tori. We consider three fiducial models of tori around rotating, both 
prograde and retrograde, and nonrotating black holes; these three fiducial models are 
also contrasted with axisymmetric simulations and a pseudo-Newtonian simulation with 
equivalent initial conditions to delineate the limitations of these approximations. There 
are both qualitative and quantitative differences in the fiducial models, with many of 
these effects attributable to the location of the marginally stable orbit, r ms (a), both 
with respect to the initial torus and in absolute terms. In the retrograde model, the 
initial inner edge of the torus is close to r ms and little angular momentum need be 
lost to drive accretion, whereas in the prograde case the gas must slowly accrete over a 
significant distance and shed considerable angular momentum. Evolution is driven by 
the magneto-rotational instability and the nonzero Maxwell stresses produced by the 
turbulence, and results in a redistribution of the specific angular momentum to near- 
Keplerian values. The magnetic energy remains subthermal within the turbulent disk, 
but dominates in the final plunging flow into the hole. The Maxwell stress remains 
nonzero in this plunging flow inside of r ms and the fluid's specific angular momentum 
continues to drop. The accretion rate into the hole is highly time variable and is de- 
termined by the rate at which gas from the turbulent disk is fed into the plunging flow 
past r ms . The retrograde model, with the largest r ms , shows the least variability in 
accretion rate. While accretion variability is a function of a, the turbulence itself is also 
intrinsically variable. A magnetized backflowing corona and an evacuated, magnetized 
funnel are features of all models. 



Subject headings: Black holes - magnetohydrodynamics - instabilities - stars: accretion 
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1. Introduction 

Numerical experiments are playing an increasingly central role in investigating accretion into 
compact objects. A particularly important application is accretion into black holes, either solar- 
mass black holes in X-ray binaries, or supermassive black holes in active galactic nuclei, quasars, and 
the like. The investigation of such systems requires the development of algorithms and simulation 
codes that solve the equations of magnetohydrodynamics (MHD) in full general relativity (GR). To 
date, however, the only global three-dimensional MHD black hole accretion simulations have been 
done using a pseudo-Newtonian potential. The pseudo-Newtonian potential is $ = —GM/{r — r g ), 
where the gravitational radius, r g , is equated with the Schwarzschild radius (Paczyhski & Wiita 
1980). This potential reproduces several important aspects of the Schwarzschild potential, in 
particular the qualitative behavior associated with the innermost stable orbit. But the system 
is still fundamentally Newtonian, and cannot account for the effects of a rotating black hole. Fully 
relativistic simulations are necessary, and in this paper we present the first such three-dimensional 
MHD black hole accretion models. 

Since the pseudo-Newtonian model is substantially simpler than full relativity, it has been the 
natural starting point for accretion simulations. We begin by briefly reviewing some of those re- 
sults. Hawley (2000; hereafter H00) computed fully global, three-dimensional prototype models of 
magnetized accretion tori using Newtonian MHD. H00 compared a variety of approximations: the 
"cylindrical disk" limit (in which the vertical component of gravity is ignored), two-dimensional 
versus three-dimensional simulations, and Newtonian versus pseudo-Newtonian potentials. This 
work demonstrated in a global context that the magnetorotational instability (MRI; Balbus &i 
Hawley 1991) drives the evolution of accretion disks by generating turbulence and transporting an- 
gular momentum through significant Maxwell stresses. Tori with near-constant angular momentum 
distributions, which served as initial conditions, are particularly unstable, and quickly evolve to 
near-Keplerian disks. It was further found that the Maxwell stress did not go to zero at the radius 
of the marginally stable orbit, r ms , the point traditionally regarded as the effective inner boundary 
of the disk. Instead, the specific angular momentum of the gas continued to decline inside r ms . 

Questions of large-scale disk structure and evolution were addressed with global MHD simu- 
lations by Hawley, Balbus & Stone (2001) and Hawley & Balbus (2002). This work examined the 
accretion of gas from a torus at 100 r g , and the subsequent formation of a disk interior to that 
point. In this model the accretion flow was assumed to be nonradiating. As a consequence, the 
resulting disk was vertically thick, with an angular momentum distribution that had a Keplerian 
slope. The flow also featured a substantial outward flowing magnetized corona. Further, the ac- 
cretion rate into the near hole region was larger than the accretion rate through the equatorial 
opening in the centrifugal barrier. Gas "backed up" at this point, creating an hot torus in the inner 
region near the hole, and an outflow along the centrifugal wall surrounding the evacuated funnel. 
The primary elements of this model — thick, nearly Keplerian disk, magnetized coronal backflow, 
hot inner torus — appear to be generic features of nonradiative accretion flows. 
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The properties of the inner region of an accretion disk are worth particularly close scrutiny 
since the most extreme relativistic effects and the greatest energy release happen near r ms . The 
properties of the inner edge were the focus of the work of Hawley & Krolik (2001, 2002) which 
used relatively high resolution simulations to study the inflow past r ms and into the black hole. 
Cylindrical disk studies investigating this region were also done by Armitage, Reynolds, & Chiang 
(2001). Krolik & Hawley (2002) examined the implications of the simulations for the inner edge of 
a black hole accretion disk. The essential points are as follows: an accretion flow must transition 
from turbulence within the disk to a plunging infall into the hole, and simulations find that the 
location of this transition is highly time- variable. Generally the transition begins outside of r ms , 
but the full transonic plunging inflow is not fully established until some point inside of T ms . The 
inner disk edge is not abrupt: average flow variables, including the stress, are relatively smooth, 
albeit time-varying, functions of radius. Magnetic stress continues to be important even within the 
plunging region. In fact, the magnetic stress becomes proportionally more important as flux-freezing 
increases the field strength relative to the gas pressure, and increases the correlation between the 
radial and azimuthal field components. 

As remarked, all the above work was done within the context of a pseudo-Newtonian potential. 
Just how well does that potential model the most important behavior in accretion flows? When is 
that approximation sufficient and when does it fail? What unique properties in the accretion flow 
are created by the spin of the hole and frame dragging, and do these produce unique observational 
signatures? Answering questions such as these requires a fully relativistic treatment. 

The development of new codes for general relativistic MHD is ongoing. Recent work includes 
that of Koide, Shibata k Kudoh (1999), Gammie, McKinney, & Toth (2003), and our own effort 
(De Villiers & Hawley 2003; hereafter DH03). In DH03 we describe the development and testing 
of the general relativistic MHD code used in this paper. The hydrodynamic portion of this code 
has already been employed (De Villiers & Hawley 2002; hereafter DH02) to study the evolution of 
accretion tori that are subject to the Papaloizou-Pringle instability (Papaloizou & Pringle 1984) 
in the full Kerr metric. Here we turn our attention to the first in a series of global simulations of 
accretion tori in the Kerr metric using the equations of MHD. Our aim with these models is to lay 
the groundwork for future, more detailed, black hole accretion simulations. 

In this paper we examine a relativistic version of the constant angular momentum torus model 
GT1 studied in H00. In §2 we write the equations of general relativistic MHD in the form that the 
code uses. We describe the constant angular momentum tori that serve as initial conditions for the 
models in this paper. We list a number of the physics diagnostics used in the simulations. In §3 we 
present the results of the three-dimensional global MHD simulations, specifically a Schwarzschild 
hole, a prograde, and a retrograde Kerr hole with near maximal spin. We contrast three-dimensional 
simulations with two-dimensional axisymmetric models. We also compare a simulation that includes 
the full range in <fi (0 < (j) < 2ir) with the results from the more economical quarter-plane simulation 
(0 < <fi < vr/2), to examine the emerging impression that quarter-plane simulations give qualitatively 
the same results, and are adequate for many studies (Hawley 2001; Nelson & Papaloizou 2003). 
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We contrast the results from the Schwarzschild simulation with an equivalent pseudo-Newtonian 
model. In §4 we review and summarize our findings. 



2. Equations of General Relativistic Magnetohydrodynamics 

We are studying the evolution of a magnetized fluid in the background spacetime of a Kerr 
(rotating) black hole. We adopt the familiar Boyer-Lindquist coordinates, (t,r,0,(f)), for which the 
line element has the form, 

ds 2 = g u dt 2 + 2 g t( f, dt d<j) + g„ dr 2 + g ee d6 2 + g^ d<p 2 . (1) 

We use the metric signature (—,+,+, +), along with geometrodynamic units where G = c = 1. 
Time and distance are in units of the black hole mass, M. The determinant of the 4-metric is g, 
and y/— g = a^fy, where a is the lapse function, a = 1/ \J —g u , and 7 is the determinant of the 
spatial 3-metric. 

The equations of general relativistic ideal MHD are the equation of baryon conservation, the 
conservation of the fluid energy-momentum tensor, and the induction equation. While the equations 
themselves are determined, there are a large number of ways that they can be expressed and solved 
numerically. In our methods paper (DH03) we derive useful forms for these equations, along with a 
set of primary and secondary variables especially suited to numerical evolution. The first of these 
is the equation of mass conservation, 

d t D + ^-d j (D^V 3 ) = 0, (2) 

where D = pW is the auxiliary density variable, p the density, and V 3 = U 3 /U 1 is the transport 
velocity, is the four- velocity, W is the relativistic gamma-factor, and U l = W/a. Spatial indices 
are indicated by roman characters i,j = 1, 2, 3. 

The internal energy equation, 

d t E + -L dj {E^yV j ) +Pd t W+^=d j (W^yV j ) = 0, (3) 

evolves E = eD, the auxiliary energy variable, where e is the specific internal energy. The isotropic 
pressure P is determined by the equation of state of an ideal gas, P = pe (V — 1), where T is the 
adiabatic exponent. For these simulations we take T = 5/3. 

The equation of momentum conservation, 

dt (S j -ab j b t ) + -^d i Sr {S 3 V*-ab ] V) + l - - a b,b)j dj + adj (^P + ^fj = 0. 

(4) 
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is written in terms of the auxiliary four-momentum, S^ = (ph + ||&|| 2 ) W U^, where h = l+e+P/p is 
the relativistic enthalpy, b^ 1 is the magnetic field four-vector in the rest-frame of the fluid, and ||6|| 2 = 
g^ v b u . The momentum is subject to the normalization condition g^ v Sn S v = —{ph + ||6|| 2 ) W 2 , 
which is algebraically equivalent to the more familiar velocity normalization = — 1 . 

Finally, the induction equation is solved in the form 

-Fa/3,7 + ^/37,a + ^7a,/3 = 0, (5) 

where F a p is the electromagnetic field strength tensor, written in terms of the constrained transport 
(CT) magnetic field variables of Evans & Hawley (1988), 

B r = F^ e , B e = F r4> ,B^ = F 0r . (6) 

We assume infinite conductivity (the flux-freezing condition), F^ v U v = 0. The induction equation 
then reads 

dj{&) =0 (i/ = 0), (7) 

8 t (B i ) - dj (V 1 B j - B i V j ) =0 {v = i), (8) 

where V^ 1 = U^/U 1 is the transport velocity with U* = W/a. The CT variables B l are related to 
the magnetic four-vector b^ by 

+ V t b t , (9) 



V47T ^jW 

and 

W 



Note that the factor a/Itt is incorporated into the definition of 6 M 



g rr V r B r + g ee V 6 B 9 + (g H V* + g t<t> ) B*\ . (10) 



For reference, the energy-momentum tensor, written using the primitive variables, has the 
familiar form 

T^ v = (ph + H&f) U"U U + ^P+M^ 9*" -b"b u , (11) 

which readily identifies the magnetic pressure, P mag = ||6|| 2 /2. 

The GRMHD code evolves time-explicit, operator-split, finite difference forms of equations 
(2) — (4) and evolves the magnetic field induction equation (8) using the CT algorithm and a mod- 
ified method of characteristic technique as described in DH03. In this paper we evolve a fiducial 
constant-angular momentum torus containing poloidal field loops orbiting around a prograde rotat- 
ing Kerr hole, a retrograde rotating Kerr hole, and a nonrotating Schwarzschild hole. The fiducial 
models use a 128 x 128 x 32 grid in (r,6,(j)). The azimuthal grid spans the quarter plane, i.e., 
< (f) < 7r/2, with periodic boundary conditions in 4>. One test simulation was run with a full 2ir 
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azimuthal grid, using 128 x 128 x 128 zones. The outer radial boundary is set to r out = 120M in all 
cases; the inner radial boundary is located just outside the horizon, with the specific value depend- 
ing on the location of the horizon as determined by the Kerr spin parameter a. The radial grid is 
spaced logarithmically to maximize the resolution near the inner boundary. The radial boundary 
conditions are zero gradient boundary conditions, where the contents of the active zones adjacent 
to the boundary are copied into neighboring ghost zones. The #-grid runs over 0.027T < 9 < 0.987T, 
with reflecting boundary conditions enforced at the polar axis. This avoids the substantial reduc- 
tion in timestep that occurs if the 9 grid runs up to the coordinate singularity at the axis. Test 
simulations have found no significant differences resulting from the reduced 9 grid. Because the 
accreting gas has angular momentum, it is excluded from the region near the rotation axis. 

A density and energy floor is employed to ensure that the dynamic range between the vacuum 
and the accreting gas does not become too great. The floor is set to be 10 -7 the initial density 
maximum in the torus. The timestep is set to one half the minimum time required for light to cross 
one grid zone. Because the grid is fixed, the timestep too is fixed throughout the simulation. 



2.1. Torus Initial State 

For the first simulations we choose a particular initial state that has been studied extensively 
in previous pseudo-Newtonian simulations: the constant-angular momentum torus containing weak 
poloidal magnetic field loops. We adopt the definition of specific angular momentum I = —U^/Ut- 
The initial axisymmetric torus is generated from the constant-/ axisymmetric GR hydrodynamic 
equilibrium equations for an adiabatic gas, as described in Hawley, Smarr, h Wilson (1984; see 
also DH02). A particular torus is determined by the choice of I, which defines a set of effective 
potential surfaces, and the specific binding energy of the effective potential corresponding to the 
torus surface, designated (t/j) . Alternately, one can specify the radius of the pressure maximum 
(where I is equal to the Keplerian value) , and the radius of the inner edge of the torus (which sets 
(Ut) s ). 

The initial magnetic field is obtained from the definition of in terms of the 4-vector 
potential, A^, 

= dp, A u - 8 u Ap. (12) 

Our initial field consists of axisymmetric poloidal field loops, laid down along isodensity surfaces 
within the torus by defining A^ = (A t , 0, 0, A^), where 



A = { k<Kf> ~ Pcut ^ f ° r 9 ~ Pcut 
* I for p < p cut 



(13) 



where p cu t is a cutoff density corresponding to a particular isodensity surface within the torus. 
Using (12), it follows that B r = —dgA^ and B e = d r A^. The constant k is set by the input 
parameter /3, the ratio of the gas pressure to the magnetic pressure, using the volume-integrated 
gas pressure divided by the volume- integrated magnetic energy density in the initial torus. In all 
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calculations reported here the average initial field strength is (3 = 100. The constant p cut is chosen 
to keep the initial magnetic field away from the outer edge of the disk. Here we use p cu t = 0.5/0 TOaa; , 
where p max is the maximum density at the center of the torus. This choice means that the initial 
field loops will be confined well inside the torus. 



2.2. Evolution Diagnostics 

One of the challenges of global, three-dimensional simulations of complex phenomena is to 
extract and distill useful diagnostic information that accurately characterizes the most important 
physical properties of the system. One computes, of course, a complete set of variables, at all grid 
zones for all timesteps; choosing from this data what to evaluate and to save is the issue. Given 
the novelty of three-dimensional MHD GR accretion simulations, there is as yet no definitive set of 
diagnostic quantities. We anticipate that standards will develop and evolve as our understanding 
of the simulation dynamics of black hole accretion flows improves. In previous pseudo-Newtonian 
calculations (e.g. H00), the evolution was characterized by several reduced quantities extracted from 
the simulation data. The set of diagnostics presented here is larger than those originally discussed 
in H00, and serves a first attempt at creating a new set of reduced quantities in a relativistic 
framework. 

Entire restart files are saved at periodic intervals, and these can be examined in detail after a 
simulation. The more difficult choice is what data to save at frequent time intervals as a running 
history of the simulation. We have used full dumps of individual variables, both fundamental and 
derived; for the present simulations we saved density, and magnetic and thermal pressure 30 times 
per orbit. We also computed a more extensive set of azimuthally-averaged variables, variables 
averaged on spherical shells, and volume-integrated quantities. We define the average of a quantity 
X on a shell at radius r as 

(X)(r) = — _ / / X^dOdcj, (14) 



A(r) 

where the area of a shell is Air) and the bounds of integration range over the 9 and (f> grids. For 
these simulations we compute shell-averaged values of density, (p), angular momentum, (pi), gas 
pressure, (P), and magnetic pressure (||fr|| 2 ). From these quantities, we derive quantities such as 
the density- weighted average specific angular momentum, (I) = (pl)/(p). 

Fluxes through the shell are computed in the same manner, but are not normalized with the 
area. We evaluate the rest mass flux (pU r ), energy flux 

(T r t ) = pU r U t + \\b\\ 2 U r U t -b r b t , (15) 

and the angular momentum flux (T r ^), saved as two separate parts, i.e., the fluid part 



(T r HF) )=pU r U 4n 



(16) 
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and the magnetic part, 

{T r H ^) = \\b\\ 2 U r U 4l -b r b^ (17) 
Again, various quantities can be subsequently derived from these fluxes and shell averages. 
Volume-integrated quantities are computed using 

[Q]=y J j Q^drdOdcj). (18) 

The volume- integrated quantities are the total rest mass, [pU*] , and total energy [T**] . In all 
the simulations these integrated values are computed and saved 30 times per orbit at the initial 
pressure maximum. 



3. Results 

3.1. Retrograde, Schwarzschild, and Prograde Black Holes 

Our fiducial simulations consider and contrast three initial torus models: a torus orbiting a 
Schwarzschild black hole, and tori orbiting rapidly rotating retrograde and prograde Kerr holes. 
Although the initial conditions cannot be made identical, we choose parameters for the initial tori 
that keep the inner edge of the disk and the location of the initial pressure maximum roughly the 
same as the black hole rotation parameter is varied. Table 1 lists the general properties of the 
models, where I denotes the specific angular momentum, rj n the inner edge of the disk (in the 
equatorial plane), (Ut) s is the binding energy of the surface of the torus, rp max the location of the 
pressure maximum (also in the equatorial plane), T or ^ the orbital period at the pressure maximum 
in units of M. For reference we also list r ms , the location of the marginally stable orbit, the values 
of {U^ms and (Ut) ms for a particle in a circular orbit at r ms (Bardeen, Press, & Teukolsky 1972), 
and T or fe( ms ), the orbital period at the marginally stable orbit in units of M. All radii in the table 
are in units of M. The retrograde SFR model is an extreme Kerr hole, a = —0.998, while the 
prograde torus orbits a rapidly rotating Kerr hole with the rotation parameter a = 0.9 chosen so 
that the marginally stable orbit is outside the static limit. 



Table 1: Global Torus Simulation Parameters. 



Model 


a 


I 


(Ut) s 


Tin 


fPmax 


Torb 


^ms 


( U 4>)ms 


( U t)ms 


T rb (ms) 


SFO 


0.000 


4.50 


-0.980 


9.5 


15.3 


376 


6.00 


3.464 


-0.942 


97.95 


SFP 


0.900 


4.30 


-0.980 


9.5 


15.4 


386 


2.32 


2.100 


-0.844 


31.04 


SFR 


0.998 


-4.80 


-0.980 


9.5 


15.8 


388 


8.99 


-4.233 


-0.962 


170.0 



Figure 1 shows the initial density profiles for these fiducial tori, plotted on the same spatial 
scale to illustrate clearly the relative size of the disks. Overlaid on the tori are contours of the initial 
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magnetic pressure. The choice of vector potential density cut-off (eq. 13), p cu t = 0.5p max , yields 
initial poloidal loops concentrated in the immediate vicinity of the pressure maximum. Although 
this keeps the field well away from the steep drop-off in density near the inner edge of the disk, the 
initial confinement of the field to the core of the disk gives rise to a violent transient as the field is 
amplified and buoyantly escapes from the core. 




Fig. 1. — Initial density and magnetic pressure profiles for the (a) SFR, (b) SFO, and (c) SFP 
models. Density is plotted on a logarithmic scale, in units of the density maximum. The colors 
cover four decades from red (high) to blue (low). Magnetic pressure, shown as overlaid contour 
lines, is also scaled logarithmically. 

The three models, SFR, SFO, and SFP, were evolved for a time equivalent to 10.0 orbits at 
their respective pressure maxima. Figures 2, 3, and 4 show polar and equatorial slices through 
the dataset for the density p taken at t = 2.0, 4.0, 10.0 orbits. During the evolution the tori 
pass through three phases: (1) the initial rapid growth of the toroidal magnetic field by shear; 
(2) the initial nonlinear saturation of the poloidal field MRI; and (3) evolution by sustained MHD 
turbulence. 

During the first phase, the shearing of the radial field quickly generates a significant toroidal 
field. The increase in magnetic pressure causes the inner region of the torus to expand, driving 
both accretion and outflow. The second phase, the violent nonlinear saturation of the MRI, begins 
at about orbit 2 and features strong radial motions. Panels (a) and (b) of each figure show plots 
at t = 2.0 orbits, where this violent transient state is at its peak. A strong, thin current sheet is 
clearly visible on the equator. A bubble of high magnetic pressure straddles the equator, driving 
lower density gas towards the black hole, while just above and below this bubble, pincers of denser 
material are also driven towards the black hole. The equatorial slices at t = 2.0 show a series of 
spirals inside r ~ 15 M, that represent a well-organized flow of material towards the hole. In model 
SFR these spiral waves reverse direction near the hole due to frame dragging. This flow feature 
was also observed in the hydrodynamic simulations of DH02. 

As gas is driven inward by magnetic pressure, it also expands upward, carrying with it signif- 
icant magnetic field. This leads to the formation of a highly magnetized outflow which moves out 
above and below the torus. This backflow is a generic feature of geometrically thick, nonradiat- 
ing accretion flows, and has been observed in previous pseudo- Newtonian simulations (Hawley & 
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Balbus 2002; Hawley & Krolik 2001, 2002). 

Through orbit 2 the flow retains a high degree of symmetry across the equator, but when 
the pincers reach the black hole, the near-horizon inflow and the inner part of the torus become 
turbulent, as can be seen at t = 4.0 orbits in panels (c) and (d). At this stage, the outer regions of 
the torus are still relatively well ordered, and two large outbound magnetic bubbles of low density 
gas can be seen above and below the equatorial plane, near the outer edge of the disk. Shortly 
after this, the entire disk becomes turbulent, and panels (e) and (f), taken at t = 10.0 orbits, are 
typical of this mature phase of the evolution. At this late stage, large spiral waves can also be seen 
moving through the outer regions of the disks. 




10 20 30 40 50 60 




30 40 50 60 
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Fig. 2. — Plots of log density for model SFR. Panels (a), (c), and (e) are polar slices through the 
disk at <j> = 0- Panels (b), (d), and (f) are equatorial slices through the disk at 6 = ir/2. Panels 
(a,b) are taken at t = 2.0 orbits; (c,d) at t = 4.0 orbits; (e,f) at t = 10.0 orbits. 
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Fig. 3. — Plots of log density for model SFO. Panels (a), (c), and (e) are polar slices through the 
disk at (j> = 0. Panels (b), (d), and (f) are equatorial slices through the disk at 8 = tt/2. Panels 
(a,b) are taken at t = 2.0 orbits; (c,d) at t = 4.0 orbits; (e,f) at t = 10.0 orbits. 
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Fig. 4. — Plots of log density for model SFP. Panels (a), (c), and (e) are polar slices through the 
disk at <j> = 0. Panels (b), (d), and (f) are equatorial slices through the disk at 6 = ir/2. Panels 
(a,b) are taken at t = 2.0 orbits; (c,d) at t = 4.0 orbits; (e,f) at t = 10.0 orbits. 

The magnetic nature of the low density bubbles at t = 2.0 and 4.0 orbits seen in figures 2-4 
is made clear by Figure 5, which shows plots of the density in the SF0 model at t = 2.0 orbits. 
Figure 5a is at t = 2.0 orbits and 5b corresponds to 4.07 orbits, a time when the bubbles are 
especially prominent. The figure shows a color plot of density with logarithmic contours overlaid 
with contours of constant magnetic pressure (dark lines). To highlight the edge of the bubbles, 
only contours of low magnetic pressure are shown (these range from 10 -4 to 10~ 3 of the maximum 
magnetic pressure on the grid). The regions inside the density-minimum bubbles correspond to 
local magnetic pressure maxima. 

The general evolutionary histories of the three models are similar: growth of a strong current 
sheet along the equator in the first orbit, initial saturation of the poloidal field MRI between the 
second and fourth orbit, and sustained turbulence thereafter. The three models also have important 
differences. The most apparent difference lies in the structure of the inflow near the horizon. At 
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Fig. 5. — Magnetic bubbles in Model SFO. Depicted are color plots of log density at (a) t = 2.0 
and (b) t = 4.07 orbits, overlaid with contours of magnetic pressure. The low density bubbles 
correspond to local maxima in magnetic pressure. 

one extreme, retrograde model SFR has a slender inflow with decreasing pressure and density, 
whereas prograde model SFP develops a new pressure maximum region inside r ~ 5 M. We refer 
to this structure as the mini-torus for its appearance when viewed in a polar slice. Model SFO is 
an intermediate case. These differences can be understood as arising from the relative distances 
between the inner edge of the initial torus and the marginally stable orbit. This distance is the 
greatest for SFR (see Table 1) and the least for SFP. 

A more detailed history of the evolution of the tori can be gained from Figure 6, a space-time 
diagram (p) (r, t) for models SFR, SFO, and SFP. The first orbit in each model shows the result of 
increasing magnetic pressure: gas is driven in toward the hole and the density in the heart of the 
torus decreases. Outgoing waves result from backflow driven from the inner region of the torus by 
this increased pressure. After orbit 2 there are strong fluctuations within the torus, and a series of 
time-variable accretion events. In models SFO and SFP there is a notable shift of the denser edge 
of the torus toward their appropriate marginally-stable orbits, and hence toward the black hole. In 
model SFR the denser material remains outside of r = 10M except for occasional accretion streams. 
After orbit 6 the evolution is less violent. The formation of the mini-torus in model SFP, which 
appears after this point in time as a high-density region inside of r ~ 5 M, is especially obvious. 

The spacetime diagram also reveals many outgoing waves. The first of these, which appear in 
the linear MRI growth stage (between and 3 orbits) originate as the backflow from the accreting 
gas. From the diagram we can measure the speed of propagation of this flow at v ~ 0.07 c. A 
prominent tilted strip, visible in each plot between t ~ 1 and 3 orbits, is associated with outward 
moving material propagating from inner edge of the inflow along the surface of the disk; this feature 
is especially prominent in animations of the data. A series of weaker outgoing waves are visible at 
later times. These are often internal spiral pressure waves rather than bulk outflows. These waves 
are also visible in animations of the data, and originate deep within the turbulent accretion flow. 
In model SFP, for example, strong waves originate within the mini-torus. 
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Fig. 6. — Space-time diagrams of shell-averaged density, (p), for models (a) SFR, (b) SFO, and (c) 
SFP. Scale is linear, and each plot is normalized to the maximum density in the initial torus. 

A quantitative view of the density distribution is provided in Figure 7, which shows (p) for 
each disk model at t = 10 orbits. In each case the density begins to decline outside of the location 
of the marginally stable orbit. In model SFP r ms = 2.32M, and the density increases roughly as 
1/r inside the initial torus's inner boundary. The region of peak density, lying between r = 3M 
and 5M, is the region referred to as the mini-torus. 




Fig. 7. — Plot of (p) as a function of radius for models SFR (dashed line), SFO (dotted line) and 
SFP (solid line) at t = 10 orbits. The density is normalized to the maximum value in the initial 
torus. 

Figure 8 shows the accretion rate through the innermost radial zones, M = —{pU r )(r m i n ,t), 
normalized to the initial torus mass for each of the models. The violent transient that results from 
the onset of the strong nonlinear phase of the MRI is clearly visible as a sharp spike at t = 2.2 
orbits in all three curves. This spike corresponds to a strong radial streaming flow of short duration. 
This type of flow is typical for the early nonlinear stage of the MRI operating on a background 
poloidal field. Strong fluctuations continue for two orbits, but by the fifth orbit all models settle 
into a steadier overall accretion rate with fluctuations about the mean on short time scales. This 
figure shows that, as a function of initial torus mass, model SFR has the greatest accretion rate, 
and model SFP the least. Accretion is easier to accomplish in model SFR because the initial torus 
inner boundary is just outside the location of the marginally stable orbit. In SFP, on the other 
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hand, the disk must work its way down almost to r = 2M before reaching the marginally stable 
orbit. This is a large enough distance that the flow must establish a new, extended turbulent disk 
structure. The mini-torus is the innermost region of that new disk. Comparing the three accretion 
rate curves also reveals that model SFR has less variability on short timescales compared to the 
other models. Variability in the accretion rate is determined by the rate at which matter is fed 
in from the turbulent disk. Fluctuations in M, therefore, originate mainly in the turbulent disk, 
outside of r ms , and the resulting time-dependence of M reflects the turbulent frequencies of the 
inner region of the disk. These frequencies are naturally lower in the SFR model. The accretion 
rates in these models were sampled only 30 times per orbit, however, which makes it difficult to 
quantify the variability, or to distinguish the variability seen in SFO from that of SFP. 




Orbits 

Fig. 8. — Accretion rate, M, at the inner radial boundary (r m j n ). The accretion rate is nor- 
malized by the initial torus mass. The value for SFR (Retrograde) is offset by 0.0006, and SFO 
(Schwarzschild) is offset by 0.0003 for clarity. 

Table 2 accounts for the total mass in each model. The values are normalized to the initial 
mass of the SFO torus. The total initial mass, Mq, and the mass at 10 orbits, Ma na i, are obtained by 
volume integration. The mass that enters the black hole, AMj n , and the mass that leaves through 
the outer boundary, AM out , are obtained by integrating the mass flux over time at the respective 
radial boundaries; the integral is approximate since the flux is sampled only 30 times per orbit. 
The fraction of initial disk mass that has left the grid through the inner boundary is given in 
the final column. Model SFR has lost the greatest percentage into the hole, and model SFP the 
least. The low mass loss through the outer grid boundary for all models means that the radial 
location of this boundary has had a minimal effect on the overall simulation. There are significant 
differences among the three models in the amount of gas moving outward, however. Model SFP 
has the greatest relative outflow; by orbit 10, 22% of the original torus mass has moved beyond 
r = 40M, a radius outside the outer edge of the initial torus. In model SFO this amount is 14%, 
and in SFR, 10%. 

The space-time diagram of gas pressure (P ga s) closely resembles that of the averaged density 
(fig. 6), and is therefore not reproduced here. The pressure decreases with time at the location 
of the original pressure maximum as the material in the torus is redistributed. The magnetic 
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Table 2: Mass Budget in Models 



Model 


M 


M fina i 


AM in 


AM out 


AM in /M 


SFR 


0.544 


0.316 


0.226 


0.0065 


0.41 


SFO 


1.000 


0.675 


0.308 


0.0246 


0.31 


SFP 


1.685 


1.323 


0.305 


0.0653 


0.18 


SFR-2D 


0.544 


0.338 


0.206 


0.0060 


0.38 


SF0-2D 


1.000 


0.709 


0.295 


0.0336 


0.30 


SFP-2D 


1.685 


1.341 


0.235 


0.1041 


0.14 


SFP-2^ 


1.685 


1.302 


0.324 


0.0661 


0.19 


PNO 


1.000 


0.637 


0.319 


0.0147 


0.32 



pressure, (P m ag), grows rapidly in the early stages of the simulations (t < 3 orbits), first primarily 
due to shear of the initial radial field, then due to the strong growth of the MRI. Once MRI-driven 
turbulence is established, the maximum of magnetic pressure shifts sharply inward with the infalling 
magnetized gas. In models SFO and SFP the location of the gas pressure maximum also moves 
toward the black hole as the turbulent disk extends inward to the radius of the marginally stable 
orbit. 

A quantitative view of the pressure evolution in the three models is given in Figure 9, where 
{Pmag) and (Pgas) are plotted as a function of time at the radius of the initial pressure maximum, 
and also near the horizon. In each case the scale is normalized to the initial pressure maximum. The 
time-evolution inside the torus is very similar in all three runs, particularly during the initial linear 
growth phase and nonlinear saturation. This indicates that the evolution of the MRI is mostly 
independent of the spin of the black hole even relatively close (r = 15M) to the hole. Gas pressure 
dominates except for a brief period around t = 2.0 orbits where the violent magnetic transient is 
at its peak. During the extended turbulent phase beyond t = 4 orbits, (3 ~ 5 to 10. Near the 
black hole magnetic pressure dominates, especially in model SFR. The two are more nearly equal 
in model SFP. 

Where the flow shifts from gas to magnetic pressure dominance is determined by the location 
of the marginally stable orbit. This is illustrated in Figure 10, which shows the averaged gas and 
magnetic pressures as a function of radius at orbit 10. Outside of r ms , the disk evolves due to 
turbulent stresses driven by the MRI, and gas pressure continues to dominate. Inside of r ms , the 
flow transitions to plunging infall, and the field evolves primarily by flux freezing. The magnetic 
pressure becomes increasingly important, helping to confine and guide the flow within the plunging 
region. The point where the averaged magnetic pressure exceeds the gas pressure is located at 
about 0.7r ms in both SFR and SFO. In SFP the inner grid boundary radius exceeds 0.7r ms , and 
the transition to magnetic dominance is not seen. 

Figure 11 shows the azimuthally averaged gas pressure, P gas , magnetic pressure, Pmag-, and 
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Fig. 9. — Panels (a), (b), and (c) show The time evolution of average pressures (P m ag) (solid line) 
and (Pgas) (dashed line) near the horizon at (a) r = 1.95 for model SFR, (b) r = 2.33 for model SF0, 
and (c) r = 1.95 for model SFP, and inside the torus at the location of the initial pressure maximum 
for (d) SFR, (e) SF0, and (f) SFP. The pressures are normalized to the initial gas pressure at the 
pressure maximum and time is in units of orbits at rp max . 



total pressure, Ptot, at i = 10.0 orbits for model SF0 (the other two models have similar profiles). 
Within the disk, total pressure is relatively smooth, consistent with a thickened equilibrium torus, 
and dominated by gas pressure. Magnetic pressure is increasingly important at the inner edge of 
the disk, and gas and magnetic pressure are comparable in the low-density coronal outflow region 
surrounding the denser disk. Magnetic pressure dominates in the funnel region near the axis, the 
pressure contours are spherically symmetric, and the magnetic field is mainly radial. The material 
in the funnel consists of fast outward-moving gas without angular momentum, several orders of 
magnitude lower in density than the maximum density in the disk (and hence below the cut-off in 
the contour plots shown above). It is important to note that the funnel region is initialized to a 
numerical vacuum, i.e. the numerical density floor of the code, and never rises much above this 
value, so conclusions about the physical properties of the funnel region must be reached cautiously. 
The point is that material from the disk, which contains significant angular momentum, is kept 
well away from the funnel by the centrifugal barrier. 

Accretion is driven by the outward transport of angular momentum through the disk via the 
MRI. Each torus began with constant specific angular momentum, /, but Figure 12 shows that by 
orbit 10 the distribution (l)(r) is approaching a more Keplerian slope. For reference, the Keplerian 
distributions of / are shown as a dotted lines. In each model the inner portion of the torus becomes 
nearly Keplerian after only one orbit of time. The outer regions evolve more slowly. The specific 
angular momentum at late time is everywhere sub-Keplerian, indicating that some radial support 
is supplied by gas pressure (e.g. fig. 11). The inflow toward the horizon shows a steady decline 
in (l) down to edge of the grid. This is indicative of continued magnetic stress acting on the fluid, 
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Fig. 10. — Transition from gas (dashed line) to magnetic pressure (solid line) inside the marginally 
stable orbit is illustrated in models SFP (top), SFO (middle), and SFR (bottom). The values are 
shell-averages at t = 10.0 orbits, and normalized to the initial torus pressure maximum value. 




Fig. 11. — Contour plots of azimuthally-averaged (a) magnetic, (b) gas, and (c) total pressure for 
model SFO (log scale). 

a result first noted in the first pseudo-Newtonian simulations (e.g., Figure 6 of H00). In the fully 
relativistic models SFR and SFO there is a noticeable roll-off in specific angular momentum for 
r < 3M (see Fig. 13). This roll-off is especially pronounced in the first orbit during the arrival 
of low angular momentum fluid that is released from the surface of the disk during a start-up 
transient. The feature remains visible and well resolved in the late stages as well. Model SFP does 
not show this abrupt drop in (I), and this is likely due to the lack of an extended plunging region 
between the marginally stable orbit and the grid inner boundary. 

A question of particular physical interest is, what is the total specific angular momentum 
carried by the accretion flow into the black hole? To answer this we divide the total angular 
momentum flux, {T r ^), by the mass flux, (pU r ), and examine the value at the inner radial boundary 
as a function of time. This is plotted in Figure 14. The values obtained by time-averages over the 
last three orbits are —4.06 for SFR, 3.32 for SFO, and 1.94 for SFP. These are all slightly below 
the value of corresponding to a circular orbit at r ms (shown as dashed lines in the figure), and 
considerably below the value in the initial torus. 

A related question is the value of the specific energy carried into the hole. This is computed 
by the ratio of the energy flux to the mass flux, {T r t) / {pU r ) at the inner radial boundary. In SFR 
this value is —0.969, in SFO it is —0.957, and in SFP it is —0.867. These are slightly less bound 
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Fig. 12. — Specific angular momentum distribution, (I) for models SFR (dot-dashed line), SFO 
(dashed line) and SFP (solid line ) at t = 10 orbits. The Keplerian distributions of angular 
momentum for all three models are shown as dotted lines. 
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Fig. 13. — Specific angular momentum distribution, (I), in the inner regions of the flows in the 
three fiducial models. The Keplerian distribution of angular momentum is shown as a dotted line; 
from top to bottom these curves correspond to the retrograde, Schwarzschild, and prograde models. 
The crosses indicate the values in the models at the radial grid zone locations. 

than the binding energy of the marginally stable circular orbit, (Ut) ms , although, in all cases, more 
bound than the value of the initial torus (—0.98). Although we do not reproduce the plot here, 
the specific binding energy of the inflow into the hole in SFR and SFO does not vary much with 
time after about 4 orbits, while the value in SFP exhibits significant fluctuations. This is another 
consequence of the absence of an established plunging inflow inside of r ms for SFP. 

The relationship between the magnetic field and the fluid's specific angular momentum is 
illustrated in Figure 15, which is an azimuthally-averaged plot of the specific angular momentum 
in the accretion flow near the black hole for model SFO. Selected contours of azimuthally-averaged 
magnetic pressure are overlaid on the plot. The specific angular momentum shows considerable 
variation with angle 9. Regions of high magnetic pressure correspond to regions of low specific 
angular momentum. The fluid with the highest specific angular momentum that reaches the black 
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Fig. 14. — Total specific angular momentum flux, (T r j,) j '(— pU r ) , into the black hole for models 
SFR, SFO, and SFP at t = 10 orbits. The value of {U^) ms for a circular orbit at r ms is shown for 
all three models dotted line. 

hole does so in a region of low magnetic pressure. The density- weighted, shell-averaged value (l) 
cannot completely capture this complexity, but it is clear from examining points at constant r near 
the inner edge of the flow, that it does fairly characterize the reduction in I that occurs near the 
horizon. 




Fig. 15. — Azimuthally-averaged contour plot of specific angular momentum {I) in the inner regions 
of the flow for model SFO. Selected contours of magnetic pressure are also shown. The magnetic 
pressure contour levels are labeled (as fraction of maximum); contour thickness is an indication of 
magnitude of ||6 2 ||. 

Much of the behavior of the three simulations depends on the location of the marginally stable 
orbit relative to the accretion torus. Krolik & Hawley (2002) defined the "turbulence edge" of a 
disk as the point where the magnetic field switches from being controlled by MHD turbulence to 
simple flux-freezing. Concomitant with this is a transition from a region where the velocity field 
is dominated by fluctuations to one where the velocity corresponds to streaming inflow. In the 
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pseudo-Newtonian simulations, Krolik & Hawley (2002) found this edge to be located at about 
1.2 r ms . Figure 16 is a density-weighted radial inflow velocity average, {—pU r )/(p), as a function 
of r/r ms at t = 10 orbits in the three fiducial models. The shell-averaged velocities within the 
turbulent portions of the disk are relatively small and variable. At a point just outside of r ms , 
however, U r begins a smooth and steady increase with decreasing radius. From r ms inward the 
slopes of the three curves show close agreement. 
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Fig. 16. — Density-weighted radial velocity average, {—pU r )/{p}, as a function of r/r ms in SFR 
(dot-dashed line), SF0 (dashed line), and SFP (solid line) at t = 10 orbits. The graph indicates 
where, in relation to the marginally stable orbit, the flow begins its transition from turbulent disk 
to streaming inflow. 



3.2. Comparisons: Axisymmetric and full 2tt Evolution 

3.2.1. Axisymmetry 

Reducing the extent of the azimuthal direction is a very effective way of reducing the compu- 
tational demands of a given problem. Axisymmetric (two-dimensional) simulations are particularly 
useful in this regard. But to what extent do axisymmetric simulations convey useful information? 
The limitations of axisymmetric pseudo-Newtonian simulations have been previously discussed by 
Hawley, et al. (2001). Here we carry out axisymmetric versions of our three fiducial models to 
investigate these issues in the general relativistic regime. These are designated SFR-2D, SF0-2D, 
and SFP-2D. The axisymmetric simulations will also provide a more complete picture of the role 
of azimuthal modes in the growth and maintenance of the MRI in the Kerr metric. 

The axisymmetric simulations are run for a time equivalent to 10 orbits at the pressure max- 
imum. As in the fiducial simulations, the magnetic pressure increases rapidly due to the growth 
of toroidal field by shear. The toroidal field MRI, however, is inherently nonaxisymmetric, so only 
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the poloidal MRI modes operate. As has previously been noted (e.g., Hawley et al. 2001), one of 
the main consequences of axisymmetry is that the evolution is dominated by a particularly violent 
form of the poloidal MRI, the so-called "channel solution" (Hawley & Balbus 1992). Figure 17 
shows the density for model SF0-2D at t = 3.0 orbits (models SFR-2D and SFP-2D are similar). 
The prominent radially-extended features are characteristic of the channel solution. The channel 
solution is itself subject to an instability which destroys its coherence, but that instability is nonax- 
isymmetric (Goodman & Xu 1994). Another known consequence of axisymmetry is that magnetic 
turbulence cannot be sustained indefinitely. In axisymmetry the toroidal field cannot be converted 
into poloidal field to create a dynamo; this is Cowling's anti-dynamo theorem. 
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Fig. 17. — Plot of log density p for model SF0-2D at t = 3.0 orbits. The axisymmetric MRI 
operating on vertical field (the "channel solution") leads to the long linear features seen here. 

The combination of the channel solution plus the anti-dynamo theorem paradoxically produces 
both more violent initial accretion and weaker long-term accretion. In the axisymmetric models the 
overall accretion totals are reduced from their three-dimensional counterparts (Table 2). However, 
in axisymmetry the initial infall of matter is more pronounced; there are significant bursts of 
accretion from the disk to the horizon. 

Figure 18 shows the accretion rate into the hole for model SFP-2D compared with that in SFP. 
The most obvious feature is that accretion in SFP-2D is highly time-variable, even episodic. In 
addition, the initial accretion flow takes longer to reach the hole; strong inflow doesn't occur until 
4.5 orbits, when there are sharp, high-amplitude spikes in M. After this, accretion continues at a 
steady low rate with periodic brief increases. In this model, matter accumulates in a mini-torus 
near the black hole, much like it does in three dimensions. At late time the mini-torus continues 
to accrete into the hole, but inflow from the main disk to the mini-torus is greatly reduced. Figure 
19 compares the averaged density (/?) at orbit 10 in model SFP and SFP-2D. One can see that the 
mini-torus is more or less distinct from the main torus, indicating that the accretion flow between 
these two regions has been greatly reduced. 

In models SFR-2D and SF0-2D (not shown) there are strong bursts of accretion between orbits 
3 and 5, and beyond that time the accretion rate declines steadily, punctuated by an occasional 
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Fig. 18. — Accretion rate, M, at the inner radial boundary (r m j n ) in models SFP-2D, SFP, and SFP- 
2-7T. The accretion rate is normalized by the initial torus mass. The value for SFP-2D (axisymmetric) 
is offset by 0.0004, and SFP (standard Prograde) is offset by 0.0002. 




Fig. 19. — Averaged density (p) as a function of radius at orbit 10 for models SFP-2D (dashed 
line), SFP (solid line) and SFP-27T (dot-dashed line). The density is normalized to the maximum 
value in the initial torus. 



flare. By the end of the three axisymmetric models, the accretion rate has quieted considerably 
and much of the matter remains out in the region of the initial torus, relatively undisturbed. 



3.2.2. Full 2vr plane 

It has been noted in pseudo-Newtonian studies that quarter-(/>-plane simulations are often 
adequate for capturing the essential aspects of accretion disk evolution (Hawley 2001; Nelson & 
Papaloizou 2003). Hawley (2001) found good qualitative agreement between the two types of grids, 
with an overall 10% reduction in field and stress levels in the restricted- simulations. 

To examine the differences between a full- and a quarter-plane simulation in the present con- 
text, we repeat SFP on a 128 3 grid with the full range in the <p coordinate (0 < <j> < 2tt). This 
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simulation is labeled SFP-2-7T. Since we simply increase the number of <j) zones by a factor of 4, the 
grid zone size, Acft, remains the same. Of the models run here, the prograde torus provided the 
greatest radial extent of turbulent disk evolution, and provides the most interesting comparison. 
Polar and equatorial slices through the dataset for the density p are similar to Figure 4, so the 
gross overall dynamics of the disk are largely unaffected by the choice of 4> grid. 

There are differences, however. A comparison of the density at orbit 10 in SFP-27T and SFP 
(fig. 19) shows that the density in the mini-torus is substantially reduced in SFP-27T. Although the 
mini-torus is initially comparable in size in the two models, it seems to partially dissipate in the late 
stages of SFP-27T. In SFP the mini-torus appears to dissipate at t « 5 orbits, but then reconstitutes 
itself. This same behavior was seen in the pseudo- Newtonian simulation of Hawley & Balbus (2002); 
an inner torus formed, dissipated, and reformed during the length of the simulation. In SFP-27T 
the density within the mini-torus declines after its formation until the end of the run. How can we 
account for this? First, the mini-torus forms, in part, because the accretion rate into the near-hole 
region doesn't necessarily match the accretion flow that passes beyond the marginally stable orbit 
and into the hole itself. The initial nonlinear MRI, for example, can rapidly dump considerable 
material into the near-hole region. At later times inflow from the extended disk may be reduced, 
while an increase in turbulence within the mini-torus can increase the rate of its accretion into the 
hole, thereby yielding the observed descrease in density within the SFP-27T mini-torus. It would 
seem, therefore, that greater stress levels in SFP-27T have a noticeable impact on disk structure. 

The long-term behavior of the accretion rate may constitute another significant difference 
between the models. Figure 18 compares the accretion rates in SFP-27T and SFP. The fluxes are 
similar in magnitude, but the full-grid simulation is smoother, and has a slightly greater mass 
flux after t = 7 orbits. Table 2 shows that the total integrated accretion in SFP-27T is about 
6% greater than that of SFP. This would indicate that the low-order azimuthal modes, whose 
existence is precluded in the quarter-grid simulations, do play a role in the accretion process and 
can increase the effective stress and angular momentum transport. The reduction in the level of the 
fluctuations, particularly when one also considers the axisymmetric result, also shows the role of 
the nonaxisymmetric modes in decreasing the coherence of the axisymmetric MRI modes (i.e., the 
channel solution). It is also possible that a reduction in A(f> would similarly reduce the fluctuations 
in M, even with a quarter-grid domain. This hypothesis remains to be tested with higher resolution 
simulations. 

3.3. Comparison with a Pseudo-Newtonian Simulation 

To compare the results from the relativistic code with those that use a pseudo-Newtonian 
potential, we ran model PNO, an analogue of SFO. These two models have several differences in 
addition to Newtonian versus relativistic physics. The PNO code is the same as used in earlier studies 
(e.g., Hawley & Krolik 2001, 2002) and, although the numerics in the two codes are based on the 
same principles of time-explicit Eulerian differencing, they are, of course, distinct algorithms. The 
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PNO code employs a cylindrical grid, (R, <fi, z), with a unit of radius twice that of the Schwarzschild 
metric, i.e., R = 2M. The radial grid has 128 zones and is equivalent to that of SFO, while the z grid 
has 128 zones, with 64 of the zones located between z = ±5. In PNO the initial torus has a pressure 
maximum at R = 7.65 and an inner edge at R = 4.75. The angular momentum is I = 3.18 (unit 
value smaller by 2 1 / 2 compared to SFO). The initial magnetic field consists of toroidal loops with 
(3 = 100, set up in the torus according to eqn. (13). The orbital period at the pressure maximum 
is T or f, = 115.6. The pseudo-Newtonian orbital period differs from that of SFO by the radial unit 
difference 2 3 / 2 (P or b oc oc r 3 / 2 ), plus a factor due to time-dilation, oc (1 — 2/r). In the GR code 
time is measured in terms of the clock of the observer at infinity, whereas there is no distinction 
between local and global time in the Newtonian computation. When comparing the results from 
the two simulations we normalize time to orbits at the initial pressure maximum. 

Figure 20 compares the history of the accretion into the hole for runs PNO and SFO. The 
plots are very similar in overall shape. In both models the initial phase of magnetic pressure-driven 
accretion, the onset of the MRI saturation phase, and the transition to more or less steady turbulent 
accretion occur at the same time, normalized by the orbital period at the pressure maximum. Table 
2 gives the mass budget for model PNO; the values are normalized to the initial torus mass. The 
mass loss through the outer boundary includes both the upper and lower z boundaries as well as 
the outer cylindrical radial boundary. PNO has lost a comparable amount of mass as the relativistic 
SFO model. 
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Fig. 20. — The accretion rate, M, at the inner radial boundary (r m j n ) in models SFO and PNO. In 
both cases the accretion rate is normalized by the initial torus mass, and time is in units of the 
orbital period at the initial torus pressure maximum. The value for PNO (pseudo-Newtonian) is 
offset by 0.0002. 

Figure 21 shows the averaged specific angular momentum (/) at 10 orbits in models SFO and 
PNO. Note that in PNO the average is taken over cylinders rather than spherical shells. However, 
since it is a density-weighted average, and most of the mass is near the equator, this difference should 
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not be too great. One sees that there is substantial agreement between the two curves, including a 
continued decline inside of the marginally stable orbit. PNO also shows a slight downturn in slope 
near the horizon, although not as pronounced as in SFO. 

5.0 
4.5 
- 4.0 

3.5 
3.0 

5 10 15 20 

r/M 

Fig. 21. — Specific angular momentum distribution, (I), for models SFO (solid line), and PNO 
(dashed line) at t = 10 orbits. The Keplerian distribution of angular momentum for SFO is shown 
as a dotted line. 

One can compute a mass averaged inflow velocity by dividing M by (p) . This averaged inflow 
velocity at 10 orbits in PNO is larger than in SFO by about 20% at r ms , and increases more rapidly 
inward. The PNO velocity reaches a maximum of about 2.4c at the horizon while SFO remains 
subluminal and approaches c. If accurate velocities are important to the analysis, a relativistic 
treatment is clearly preferred. 

4. Discussion 

In this paper, we present a preliminary survey of the properties of accreting tori in the Kerr 
metric. This work is the relativistic analogue of earlier work by Hawley (2000) which used a 
pseudo-Newtonian potential. We carried out three-dimensional simulations of an accretion torus 
around prograde and retrograde Kerr holes with large a values, as well as around a nonrotating 
Schwarzschild hole. These simulations are contrasted with axisymmetric simulations and a pseudo- 
Newtonian simulation with equivalent initial conditions. 

In comparing the three fiducial models, we find that the most important effect of black hole 
rotation is in setting the location of the marginally stable orbit. In retrograde model SFR, which 
has the outermost r ms , the initial torus begins with its inner edge almost at the marginally stable 
orbit. As pressure builds and angular momentum is transported out through the torus, a slender 
accretion flow is set up which remains relatively smooth in the late stages of the simulation. The 
ease with which the matter is accreted is indicated by the observation that over the course of the 
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simulation, 42 % of the initial disk mass enters the black hole in 10 orbits. At the other extreme, in 
the prograde model, SFP, the gas must slowly accrete from the initial torus down to r ms = 2.32 M, 
a radius just outside the static limit. This accretion is driven by angular momentum transport 
from the MHD turbulence, and the flow develops into a dense, slightly sub-Keplerian disk. The 
inner edge of this disk features a local pressure maximum inside r w 5 M which we refer to as the 
"mini-torus." All of these factors reduce the total accretion through the inner radial boundary into 
the hole over the course of the 10 orbit evolution. Only 18% of the initial torus mass is lost into 
the hole. A greater fraction of the mass is lifted to radii outside the outer boundary of the initial 
torus. The Schwarzschild model SFO is an intermediate case which develops both a thick, turbulent 
accretion flow and a plunging inflow. SFO loses 31 % of its initial mass to the black hole. 

In all three models the accretion rate into the hole is highly time variable. SFR shows less 
variability on short time scales. The accretion rate into the hole is determined by the rate at 
which material from the turbulent disk is fed into the plunging flow from the disk's "turbulence 
edge," a point located near, but generally slightly outside, the radius of the marginally stable 
orbit. The accretion rate variability will be roughly determined, therefore, by the orbital frequency 
at the turbulence edge. Assuming that fluctuations in M have observable consequences (e.g., as a 
possible origin of high-frequency QPOs), then these simulations confirm the idea that the hole's spin 
parameter a helps determine the variability frequency. The intrinsic variability of the turbulence 
itself, however, prevents the relationship between a and variability frequency from ever being more 
than approximate. 

As in our earlier hydrodynamic simulations, we observe the effect of frame dragging most 
prominently in the retrograde model SFR, where the spiral accretion stream shows the characteristic 
"dog leg" pattern near the hole as the stream is dragged into co-rotation with the black hole inside 
r rs 3 M. The effects of prograde frame dragging are much more difficult to discern directly. 

We find that, as in earlier pseudo-Newtonian simulations, an extended magnetized corona 
forms where the magnetic pressure is comparable to the gas pressure. Much of this corona is 
flowing radially outward. The largest outflows occur in SFP; again, because it has higher specific 
angular momentum relative to the marginally stable orbit, the gas is not as easily lost to the hole 
as in SFR. 

In all the models the initial magnetic field corresponded to a volume-averaged strength of 
(3 = 100. In the subsequent evolution the disk remains dominated by gas pressure, with (3 ~ 10, 
while the coronal regions are more near equipartition, with (3 ~ 1. Magnetic pressure dominates 
within the centrifugal funnel region near the black hole. However, there is no significant matter 
density within the funnel; the specific angular momentum of the gas prevents the accreted gas from 
approaching the axis. 

An observation that clearly emerges from the simulations is the dominant role of the magnetic 
fields within the plunging inflow inside the marginally stable orbit. In models SFO and SFR there 
is a noticeable transition from a gas-pressure dominated to a magnetic-pressure dominated flow 
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inside r « 0.7 r ms . For model SFP, however, this transition does not occur; the marginally stable 
orbit lies too close to the inner radial boundary. Instead we see near-equality of the two pressures 
at the inner radial boundary. A higher-resolution simulation of model SFP which includes more 
zones within the ergosphere is required to resolve the plunging inflow. 

The magnetic fields within the plunging region continue to exert a stress on the fluid. Regions 
of strong field correspond to the lowest specific I values in the gas entering the hole. In all three 
models the specific angular momentum flux entering the hole, (T r < p)/M, is slightly less than the 
value associated with the marginally stable orbit: 96% for SFO and SFR, and 92 % for SFP. 

We also compared the fiducial simulations to axisymmetric simulations. The overall initial 
evolution of the axisymmetric and three-dimensional simulations are similar. The magnetic pres- 
sure increases, causing the torus to expand and creating outflows, while the angular momentum 
within the torus is redistributed to a more Keplerian form. However, the duration of MRI-driven 
turbulence is short-lived in the axisymmetric case since the azimuthal modes that sustain the MRI 
are absent under axisymmetry. Further, as the SFP-2D simulation demonstrates, the turbulence 
can die out in one region of the flow while continuing in others, therefore significantly altering the 
global evolutionary properties. 

As noted earlier by Hawley & Balbus (2002), axisymmetric simulations tend to be dominated 
by the poloidal field MRI channel solution which is characterized by radially extended filamen- 
tary structure and large angular momentum transport. The channel solution is also present in 
the three-dimensional simulations, but it is short-lived due to nonaxisymmetric modes, and its 
breakup heralds the onset of MRI-driven turbulence. These axisymmetric MRI modes are highly 
time- variable, causing accretion to occur in bursts. Comparing the results of the full 2ir plane simu- 
lation SFP-27T with the standard quarter-plane simulation shows that with greater axial extent, the 
accretion flow into the hole becomes smoother in time. This points to the role of nonaxisymmetric 
modes in moderating the influence of the axisymmetric channel solution. 

Since axisymmetric simulations cannot model essential aspects of the accretion process, their 
utility is somewhat limited. Used with care they can provide an inexpensive way to evolve an 
idealized initial condition, such as a constant-! torus with loops of field, into a more extended near- 
equilibrium structure that then serves as the starting point for a full three-dimensional simulation. 
They also provide a valuable tool for surveys of various initial conditions. Since axisymmetric 
simulations can be much more highly resolved, they also should continue to be useful as numerical 
experiments to investigate specific physical processes in detail. 

In comparing full- and quarter-plane simulations we find that, while there are noticeable differ- 
ences, the quarter-plane simulations capture much of the qualitative behavior. This is the consistent 
with the conclusion from previous pseudo-Newtonian studies. Their affordability makes them very 
useful for running a range of models, such as the present work. The accretion rate in the full-plane 
simulation is less strongly variable, presumably due to a reduction in the power of the axisymmetric 
poloidal field MRI. The time-averaged accretion rate, however, is higher, and this influences the size 
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of the mini-torus that forms in the prograde runs. The optimal strategy may be to run models on 
restricted grids until the startup transients are past and the turbulent disk is well-established. The 
0-grid can then be increased to get an improved representation of the approximately steady-state 
disk. 

A direct comparison of relativistic model SFO and pseudo-Newtonian model PNO reveals good 
agreement: the two simulations show similar late-time angular momentum distributions as well 
as comparable accretion histories through the inner radial boundary. In the regions outside the 
marginally stable orbit relativistic and pseudo-Newtonian results overlap, which should not be 
surprising. This is a region of weaker gravitational field and the pseudo-Newtonian potential is a 
very good approximation to the Schwarzschild metric. The differences between the two simulations 
arise in the region inside r ~ 3M, where general relativistic effects become prominent. This is also 
the region where the velocities in pseudo-Newtonian models are likely to be the most inaccurate. In 
relativity velocities are limited to less than c, whereas there is no such limit in a pseudo-Newtonian 
simulation. We further note that there is no unique time-normalization that is appropriate for all 
locations in the flow when comparing a relativistic simulation with a Newtonian one. If details of 
the time-variability of the flow in the near hole region are important, relativistic calculations are 
required. Simulating accretion into Kerr holes also requires a relativistic treatment. 

As the time-dependent models improve, it is important to connect the results to observable 
phenomena. Recently, Armitage & Reynolds (2003) took a first step toward this goal with a study 
of the effects of time variability on certain emission processes in the inner regions of a simulated 
black hole accretion disk. Machida & Matsumoto (2003) have evolved a global disk to search for 
regions within the plunging flow that might account for X-ray flares. These studies must ultimately 
be done with full relativity. In the present work we see both qualitative and quantitative differences 
among the models with different a values in the plunging region inside the marginally stable orbit, 
where effects of black hole rotation are increasingly important. More detailed analyses of this 
plunging region will be the subject of subsequent papers. 

This work was supported by NSF grant AST-0070979 and NASA grant NAG5-9266. The 
simulations were carried out on the Origin 2000 system at NCSA, and the Bluehorizon system of 
NPACI. 
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